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This paper is considered with joint estimation of state and time-varying noise covariance matrices in non-linear stochastic state 
space models. We present a variational Bayes and Gaussian filtering based algorithm for efficient computation of the approximate 
filtering posterior distributions. The Gaussian filtering based formulation of the non-linear state space model computation allows 
usage of efficient Gaussian integration methods such as unscented transform, cubature integration and Gauss-Hermite integration 
along with the classical Taylor series approximations. The performance of the algorithm is illustrated in a simulated application. 
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1. Introduction 

In this paper, we consider Bayesian inference on the state x k 
and noise covariances Y, k in heteroscedastic non-linear stochas- 
tic state space models of the form 



x*~N(/(**_i), Q k ) 
y k ~N(h(x k ),-L k ) 
~ p(L k |Zi-i), 



(1) 
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where x k e W is the state at time step k, and y k e M. is the 
measurement, Q k is the known process noise covariance and Ej 
is the measurement noise covariance. The non-linear functions 
/(•) and h(-) form the dynamic and measurement models, re- 
spectively, and the last equation defines the Markovian dynamic 
model for the dynamics of the unknown noise covariances Ej. 

The purpose is to estimate the joint posterior (filtering) dis- 
tribution of the states and noise covariances: 



P(x k ,^ k \yi: k ), 



(2) 



where we have introduced the notation y\± — yi, . . . ,y k . 

If the parameters Q k and Y. k in the model Q were known, 
the state estimation problem would reduce to the classical non- 
linear (Gaussian) filtering problem 0] [2]]. However, here we 
consider the case, where the noise covariances Y. k are unknown. 
The formal Bayesian filtering solution for general probabilistic 
state space models, including the one considered here, is well 
known (see, e.g., 0]]) and consist of the Chapman-Kolmogorov 
equation on the prediction step and Bayes' rule on the up- 
date step. However, the formal solution is computationally in- 
tractable and we can only approximate it. 

In a recent article, Sarkka and Nummenmaa [3| introduced 
the variational Bayesian adaptive Kalman filter (VB-AKF), 
which can be used for estimating the measurement noise vari- 
ances along with the state in linear state space models. In this 
paper, we extend the method to allow estimation of the full 



noise covariance matrix and non-linear state space models. The 
idea is similar to what was recently used by Piche et al. (4) in 
the context of oulier-robust filtering, which in turn is based on 
the linear results of J3] . 

We use the Bayesian approach and use the free form varia- 
tional Bayesian (VB) approximation (see, e.g., 0|7][8 1) for the 
joint filtering posteriors of states and covariances, and the Gaus- 
sian filtering approach ESj for handling non-linear models. 
The Gaussian filtering approach allows us also to utilize more 
general methods such as unscented Kalman filter (UKF) ifTTl . 
Gauss-Hermite Kalman filter (GHKF) [9 1, cubature Kalman fil- 
ter (CKF) 0~2), and various others Ifl3l [T4l [T5l along with the 
classical methods fl~l|2|. 

The variational Bayesian approach has been applied to pa- 
rameter identification in state space models also in lfT6l [T71 [181 
and other Bayesian approaches are, for example, Monte Carlo 
methods Ifl9l l20l ETIl and multiple model methods |) 22ll . It is 
also possible to do adaptive filtering by simply augmenting the 
noise parameters as state components |2| and use, for example, 
above-mentioned Gaussian filters for estimation of the state and 
parameters. 

1.1. Gaussian Filtering 

If the covariances in the model ([TJ were known, the filtering 
problem would reduce to the classical non-linear (Gaussian) op- 
timal filtering problem Q]|2|. This non-linear filtering problem 
can be solved in various ways, but one quite general approach 
is the Gaussian filtering approach [2 9 10], where the idea is to 
assume that the filtering distribution is approximately Gaussian. 
That is, we assume that there exist means m k and covariances 
P k such that 

p(x k | yi-k) ~ N(x k I m k , P k ). (3) 

The Gaussian filter prediction and update steps can be written 
as follows |9|: 
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Prediction: 



2. Main Results 



J (/(in) - m£) (/(in) - m- k f (4) 



fc— i 



(5) 



Update: 

jU* = J" N(xa; | m^, dx t 

5* = J(h(x k ) - fi k ) (h(x k ) - p k ) T 

xN(x k \m k ,P k )dx k + I, k 

Ck = J (x k - m~) (h(x k ) - p k ) T N(xk I m k , P k ) dx k 

K k -C k S k x 

m k = m k +K k (y k - fi k ) 

P k =P k -K k S k K T k . 



With different selections for the Gaussian integral approxima- 
tions, we get different filtering algorithms IflOl . 

1.2. Variational Approximation 

In this paper, we approximate the joint filtering distribution 
of the state and covariance matrix with the free-form variational 
Bayesian (VB) approximation (see, e.g., ll6ll7l[8l [T6l '): 



p{x k ,^k \yv.k) ~ Qx(Xk) Qx(Zk), 



(6) 



where Q x (x k ) and Qz(I. k ) are the yet unknown approximating 
densities. The VB approximation can be formed by minimizing 
the Kullback-Leibler (KL) divergence between the true distri- 
bution and the approximation: 

KL[<2*(*jO Qsffik) II p(x k , -L k I yi:*)] 



Qx(x k ) (h&k) log — — — r dx k dSfc. 

p(Xk, U I y\:k) 



Minimizing the KL divergence with respect to the probability 
densities, we get the following equations: 



■If 



Qx{x k ) oc exp logp(y k ,x k ,'£k\yi:k-i)Qz(?'k)$2'k 



a 



Qx<£k) K exp \ogp(y k ,x k ,l. k \yi: k -i)Q x {x k )dx k 



(1) 



These equations can be interpreted and used as fixed-point iter- 
ation for the sufficient statistics of the approximating densities. 

In the original VB-AKF [3], the VB approximation was de- 
rived for linear state space models with diagonal noise covari- 
ance matrix. In this paper, we generalize it to non-linear sys- 
tems with non-diagonal noise covariance matrix. 



2.1. Estimation of Full Covariance in Linear Case 

We start by considering the linear state space model with un- 
known covariance as follows: 



p(x k \x k -i) = N(x k \A k Xk-i,Qk) 
p(y k | x k , E*) = N fe I H k x k , 



(8) 



where A k and H k are some known matrices. We assume that the 
dynamic model for the covariance is independent of the state 
and of the Markovian form p(T. k an d set some restric- 

tions to it shortly. In this section we follow the derivation in 
0, and extend the scalar variance case to the full covariance 
case. 

Assume that the filtering distribution of the time step k - 1 
can be approximated as product of Gaussian distribution and 
inverse Wishart (IW) distribution as follows: 

p(Xk-u^k-\ \ yv.k-i) = 

Nfrk-^mk-uPk-dWVk-ilvk-uVk-!). 

where the densities, up to non-essential normalization terms, 
can be written as 11231 : 

1 

IW(E|v,V)cc |2r (v+M+1)/2 exp|-^tr(yi- 1 )J. 

That is, in the VB approximation (rol, Q x (x k ) is the Gaussian 
distribution and Qz(L k ) is the inverse Wishart distribution. 

We now assume that the dynamic model for the covariance 
is of such form that it maps an inverse Wishart distribution at 
the previous step into inverse Wishart distribution at the current 
step. That is, the Chapman-Kolmogorov equation [1| gives the 
prediction 



N(x \m,P)oc \P\- 112 exp \-^{x - mf P l (x - m) | 



p(^k\yi-.k-\) 

= j p^lS^IWffinlviM.Vi-OdZn 
= IW(Z,|v^,Vp, 



for some parameters v, and V, . We postpone the discussion on 



how to actually calculate these parameters to Section 2.2 For 
the state part we obtain the prediction 



P(Xk\yi:k-l) 

N(x k \A k x k -u 8t)N(in \ m k -i,P k -\)&x k -\ 
N(x k \m k ,P k ), 



where m, and P. are given by the standard Kalman filter pre- 
diction equations: 



: A k m k -i 
A k Pk-iA T k +Q k 



(9) 



Because the distribution and the previous step is separable, and 
the dynamic models are independent we thus get the following 
joint predicted distribution: 

p(x k ,Y k \y v . k - 1 ) = 

N(x k \m k ,P k )lW(Z k \v k ,V k ). 

We are now ready to form the actual VB approximation to the 
posterior. The integrals in the exponentials of (|7]i can now be 
expanded as follows (cf. 01): 



log p(y k , x k , | yi;jfc_i) dL k 



= ~~(yk - H k x k ) T {\ l )i(y k - H k x k ) 



I 



(10) 



- ^ (** - m k ) T (P k ) (x k -m k ) + C\ 
log p(y k , **, E t | y\±-\) Q x (x k ) dx k 
= -I(v- + „ + 2) logis.i-itrfy;!, 1 } 



:((j>k ~ H k x k ) T *L k x (y k - H k x k )) x + C 2 , 



where <-) s = / (•) &(£,) 6Z k , <■>, = / (•) Q x (x k ) dx k , and C u C 2 
are some constants. If we have that Qz(L k ) = IW(Lt [ v k , V k ), 
then the expectation in the first equation of (fTOb is 



&th = (yk-n-i)vz\ 



(ID 



Furthermore, if Q x (x k ) = N(x k [ m k , P k ), then the expectation in 
the second equation of ( |T0] > becomes 

<(y* - H k x k ) T YT k 1 Oi - //i 

= tr{^P,HjS^} (12) 
+ tr {(>>£ - mi) (y* - m k ) T lr k 1 } . 

By substituting the expectations ( fTT| and (jT2]» into ( p~0] > and 
matching terms in left and right hand sides of |7]i results in the 
following coupled set of equations: 

S k = H k p- k H T k +{v k -n-\r l V k 
K k = P k H T k S k 
m k = m k + K k (y k - H k m~ k ) 
P k =p- k -K k S k K T k 
v k = v k + 1 

V k = V k + H k P k H T k + (y k - H k m k ) (y k - H k m k ) T . 

The first four of the equations have been written into such sug- 
gestive form that they can easily be recognized to be the Kalman 
filter update step equations with measurement noise covariance 
(n-n-1)" 1 V k . 

2.2. Dynamic Model for Covariance 

In analogous manner to |3), the dynamic model p(T, k \ Sh) 
needs to be chosen such that when it is applied to an inverse 



(13) 



Wishart distribution, it produces another inverse Wishart distri- 
bution. Although, the explicitly construction of the density is 
hard, all we need to do is to postulate a transformation rule for 
the sufficient statistics of the inverse Wishart distributions at the 
prediction step. Using similar heuristics as in 0, we arrive at 
the following dynamic model: 



v k =p(v*_i -n 



!) + « + ! 



(14) 



where p is a real number < p < 1 and B is a matrix 
< \B\ < 1. A reasonable choice for the matrix is B — -\fpl,in 
which case parameter p controls the assumed dynamics: value 
p = 1 corresponds to stationary covariance and lower values 
allow for higher time-fluctuations. The resulting multidimen- 
sional variational Bayesian adaptive Kalman filter (VB-AKF) 
is shown in Algorithm[T] 



Predict: Compute the parameters of the predicted distribu- 
tion as follows: 

m k = A k m k -i 
P k =A k P k _ 1 A T k +Q k 
v k =p (y k -\ - «-l)+n + l 
V7 = B V k -x B T , 



Update: First set m 



and V 



(0) 



- V k and the iterate the following, say N, steps 
i = l,...,N: 



(0) 



1 + vl 



Sf l) =H k P7H T k+ {v k 



i)- 1 vf 



K 



<i+i) 



= p- k Hi[sr ) r l 



in,. 



■tf M) (y k - 



H k m k ) 



k 

' k ^k 

' [) = V k + H k Pf H[ + (y k - H k mf) (y k - H k mff 



pf l) = Pk-K 



V, 



Algorithm 1: The multidimensional Variational Bayesian 
Adaptive Kalman Filter (VB-AKF) algorithm 



2.3. Extension to Non-Linear Models 

In this section we extend the results in the previous section 
into non-linear models of the form ([T}. We start with the as- 
sumption that the filtering distribution is approximately product 
of a Gaussian term and inverse Wishart (IW) term: 

P(**-I,2jfc_i lyiA-i) = 

NO^Im^i.P^OIWCXi.iln-i, V w ). 

The prediction step can be handled in similar manner as in the 
linear case, except that the computation of the mean and co- 
variance of the state should be done with the Gaussian filter 



3 



prediction equations Q instead of the Kalman filter prediction 
equations i9). The inverse Wishart part of the prediction re- 
mains intact. 

After the prediction step, the approximation is then 

p(jCjfe,Sit|yi:t_i) = 

N(x k \m- k ,p- k )W(Z k \v k ,V k ). 

The expressions corresponding to ( p~0] > now become: 



log p(y k , Xk, 2* |yi : *-i) Ge(S0 dZi 

= -^Cy* - A(x*)) r <SJ 1 > s (y* - *(**)) 

- ^(** - i%f (P^j 1 (x* - mT k ) + d 
log pCy*;, jc t , | yi : n) Q x (x k ) dx k 
= -I(y- + „ + 2) logEkl-itrfv^^ 1 } 



(15) 



2 
1 

2* 



- ^<(y* - Kx k )) T ^ (y k - h(x k ))) x + C 2 . 



The expectation in the first equation is still given by the equa- 
tion ( fTTj ), but the resulting distribution in terms of x k is in- 
tractable in closed form due to the non-linearity h(x k ). For- 
tunately, the approximation problem is exactly the same as en- 
countered in the update step of Gaussian filter and thus we can 
directly use the equations |5]l for computing Gaussian approxi- 
mation to the distribution. 

The simplification ( p"2| ) does not work in the non-linear case, 
but we can rewrite the expectation as 



{(y k -h{x k )) T Y k l (y k -h{x k ))) x 



= tr (<(y t - h(x k )) (y k - h(x k )f) x Z k 1 }, 



(16) 



where the expectation can be separately computed using some 
of the Gaussian integration methods in 1 10 1 . Because the result 
of the integration is just a constant matrix, we can now substi- 
tute (JTTJ and ( fT6| i into ([15) and match the terms in equations |7) 
in the same manner as in linear case to obtain the equations: 

Uk = J Kx k ) N(x k | m\~, Pp dx k 

S k = J (Kxk) ~ Uk) (h(Xk) ~ Hkf 

X N(x k | m k ,P k ) dx k + (v k -n-\y l V k 
C k = J (x k - m~) (h(x k ) - /i k ) T N(x k | m k , P k ) dx k 



K k - C k S k 

m k = m k + K k (y k - Uk) 
P k = P k -K k S k K T k 
v k = v~ k + 1 

V k = V k - + J (y k - h(x k )) (y k - h(x k )) 7 
xN(x k \m k ,P k )dx k . 



(17) 



2.4. The Adaptive Filtering Algorithm 

The general filtering method for the full covariance and non- 
linear state space model is shown in Algorithm[2] Various use- 
ful special cases and extensions can be deduced from the equa- 
tions: 

• The Gaussian integration method will result in differ- 
ent variants of the algorithm. For example, the Taylor 
series based approximation could be called VB-AEKF, 
unscented transform based method VB-AUKF, cubature 
based VB-ACKF, Gauss-Hermite based VB-AGHKF and 
so on. For the details of the different Gaussian integration 
methods, see, HOI HU HI El El [11113. 

• The linear special case can be easily deduced by compar- 
ing the equations ( [T3) and ( fP7| > to the equations in the algo- 
rithm. That is, one simply replaces the Gaussian integrals 
with their closed form solutions. 

• The diagonal covariance case, which was considered in 
0, can be recovered by updating only the diagonal ele- 
ments in the last equation of the Algorithm and keeping 
all other elements in the matrices V ( k zero. Of course, 
the matrix B in the prediction step then needs to be diag- 
onal also. Note that the inverse Wishart parameterization 
does not reduce to the inverse Gamma parameterization, 
but still the formulations are equivalent. 

• Non-additive dynamic models can be handled by simply 
replacing the state prediction with the non-additive coun- 
terpart. 



3. Numerical Results 

3.1. Range-Only Tracking in a Non-homogeneous Noise Field 

In this simple example we illustrate the performance of the 
developed adaptive filters by tracking a moving target with 
sensors, which measure the distances to the target moving in 
2-dimensional (u, v) space. The measurements are corrupted 
with noise having time-varying correlations between the sen- 
sors. The correlations arise, because the noise in the measure- 
ments is caused by localized variations in the environment and 
when the spatial paths of the measured signals are similar, the 
noises are correlated. 

The state is contains the position and velocity of the target 
x k = (u k v k u k v k ) T and the dynamics of the target are modeled 
by the standard Wiener velocity model. The distance measure- 
ments from m sensors read 



y k 



^(s l u -u k ) 2 + (s[, 



v k ) 2 + r' i=l, 



(18) 



where (s' u , s[,) is the position of ith sensor and ri is the ith com- 
ponent of a Gaussian distributed noise vector r k ~ N(Q, Z^). In 
this experiment the noise to each distance measurement is gen- 
erated by drawing a random sample from a discretized Gaussian 
random field z k e W z and then collecting all the values of the 
field connected to the line between the sensor and the target. 
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We take the time-white continuous-valued random field z(u, v) 
to be zero mean and to have the covariance function 



k(u, v, u' , v') = cr\ g S(u — u')6(v — v') 



y| l,4 ; (w, v) v') £,(«, v, w', v'), 



(19) 



where the first term corresponds to white background noise 
and the latter terms to additive correlations for points inside 
bounded regions A, c K 2 , i - 1 . . . , «a with covariance func- 
tions kt(u, v, u', v'). We set the covariance functions to 



ki(u, v, u , v') — cr? 1 5(u - u) 6(v - v') 



+ °"ma R n,, eX P 



( 1 N 

j^u-u'f + (v-v') 2 ) 



(20) 



which means that noise inside the bounded regions consists of 
independent (white) and correlated components. 

The simulation scenario is illustrated in Figure [T] For the 
lightly shaded area (Aj) the covariance function parameters 
were cr 2 . = 0.01 2 , cr 2 im ,. = 0.1 2 and /; = 2, and inside the 
darkly shaded area (A 2 ) . = 0.01 2 , o^ agn ,. = 0.2 2 and 1, = 2. 
The variance of the background noise was set to cr 2 g = 0.01 2 . 
The spectral density of the process noise was set to q — 2 and 
the time step to T — 0.01. The trajectory shown in Figure [T] 
was discretized to 1000 time steps and then measurements were 
generated according to the procedure described above. Given 
the measurements, the target was tracked with the following 
methods: 

• UKF-t: Unscented Kalman filter with measurement co- 
variance set to true value on each time step. 

• UKF-o: Unscented Kalman filter with fixed diagonal mea- 
surement covariance matrix with different standard devia- 
tions o~ = 0.1, . . . ,3. 

• VB-AUKF-f: The proposed adaptive filter with ADF ap- 
proximations made with UKF. The parameter p in dynamic 
model of measurement noise was set to p = 1 - exp(-3). 

• VB-AUKF-d: Same as VB-AUKF-f with the exception 
that the measurement covariance is forced to be diagonal. 

With all methods the parameters of UKF was set to default val- 
ues a = = 0,/f = 3 — m. The RMSE values in tracking 
the position of the target with the tested methods are shown 
in Figure [2] for a typical simulation run. It can be seen that the 
best results can be achieved with exact measurement covariance 
while the estimation of full covariance improves the results of 
VB-AUKF over the diagonal case. Obviously, UKF with fixed 
diagonal measurement covariance is clearly the worst of the all 
the tested methods. Figure[3]shows the estimates of the element 
of E£ produced by VB-AUKF-f together with their true values. 
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Figure 1: Simulation scenario in the Range-Only Tracking Demonstration. Cir- 
cle denotes the starting location of the target, triangles the locations of the sen- 
sors and the dashed line the trajectory of the target. Inside the shaded areas the 
noise field has spatial correlations. 




Figure 2: Range-Only Tracking: RMSE values for tracking the position of the 
target on a typical simulation run. The different standard deviation values used 
for UKF-o are on X-axis. 



3.2. Multi-Sensor Bearings Only Tracking 

As an example, we consider the classical multi-sensor bear- 
ings only tracking problem with coordinated turning model, 
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estimates can be obtained by tuning p to be larger, but that introduces more lag x ^ g 

to estimates. 

Figure 4: The simulated trajectory and the estimate obtained with VB-ACKF. 



where the state x = (u, u, v, v, a>) contains the 2d location (u, v) 
and the corresponding velocities (u, v) as well as the turning 
rate to of the target. The dynamic model and the measurement 
models for sensors i = 1, . . . ,4 are given as: 



yk 



1 


sin(tL»A/) 








cos(wAf) 








l-cos(wAz) 


1 





sin(wAf) 





lo 










/ v k - si 




arctan r 


)♦ 




\uk - 4 





_ ^ 1 -cos(q)AQ ^ 

- sin(dLtAf) 

sin(o»A/) 

cos(wAf) 




0^ 




U 



x k-\ + 1k-\ 



(21) 



1. 



where ~ N(0, Q) is the process noise and = (r^i, . . . , r^) 
are the measurement noises of sensors with joint distribution 

~ N(0, T^k), where S^. is unknown and time varying. 

We simulated a trajectory and measurements from the model 
and applied different filters to it. We tested various Gaussian 
integration based methods (VB-AEKF, VB-AUKF, VB-ACKF, 
VB-AGHKF) and because the results were quite much the same 
with different Gaussian integration methods (though VB-AEKF 
was a bit worse than the others), we only present the results ob- 
tained with VB-ACKF. Figure [4] shows the simulated trajectory 
and the VB-ACKF results with the full covariance estimation. 
In the simulation, the variances of the measurement noises as 
well as the cross-correlations varied smoothly over time. The 
simulated measurements are shown in Figure [5] 

Figure [6] shows the root mean squared errors (RMSEs) for 
the following methods: 

• CKF-t: CKF with the true covariance matrix. 

• CKF-o: CKF with a diagonal covariance matrix with di- 
agonal elements given by the value on the x-axis. 

• VBCKF-f. CKF with full covariance estimation. 

• VBCKF-d: CKF with diagonal covariance estimation. 




200 400 600 800 1000 

Figure 5: The simulated measurements. 

As can be seen from the figure, the results of filters with co- 
variance estimation are indeed better than the results of any fil- 
ter with fixed diagonal covariance matrix. The filter with the 
known covariance matrix gives the lowest error, as would be 
expected, and the filter with full covariance estimation gives a 
lower error than the filter with diagonal covariance estimation. 

4. Conclusion and Discussion 

In this paper, we have presented a variational Bayes and 
Gaussian filtering based algorithm for joint estimation of state 
and time-varying noise covariances in non-linear state space 
models. The performance of the method has been illustrated 
in simulated applications. 

There are several extensions that could be considered as well. 
For instance, we could try to estimate the process noise covari- 
ance in the model. However, it is not as easy as it sounds, be- 
cause the process noise covariance does not appear in the equa- 
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Figure 6: Root mean squared errors (RMSEs) for different methods. 

tions in such simple conjugate form as the measurement noise 
covariance. Another natural extension would be the case of 
smoothing (cf. |4|). Unfortunately the current dynamic model 
makes things challening, because we do not know the actual 
transition density at all. This makes the implementation of a 
Rauch-Tung-Striebel type of smoother impossible — although 
a simple smoothing estimate for the state can be obtained by 
simply running the RTS smoother over the state estimates while 
ignoring the noise covariance estimates completely. However, 
it would be possible to construct an approximate two-filter 
smoother for the full state space model, but even in that case 
we need to put some more constraints to the model, for exam- 
ple, assume that the covariance dynamics are time-reversible. 
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• Predict: Compute the parameters of the predicted distribu- 
tion as follows: 

ml = J/(l(-l)N(ln|ffln,Pn)dln 

P k=J (f( x k-i) ~ m k ) (/(**- 1) - m k ) T 

xN(in \ m k -i,Pk-i)dx k -i + Q k 
v k = p{v k -\ - n-\)+n+\ 
V k =BV k - X B T , 

• Update: First set rn^ = m k , f[ 0) = P k , v k = 1 + v k , and 
V^ 0) = V k and precompute the following: 

Hk = J li(x k )N(x k \m k ,P k )dx k 

(h(x k ) - Hk) (Kxk) - Hk? 
x N(x k \m k ,P' k )dx k 
C k = J(x k - m~) (h(x k ) - Hk) T 
xN(x k \m k ,P k )dx k 
Iterate the following, say N, steps i = l,...,N: 

Sf l) = T k + {y k -n-\T l vf 

4 +1) =C t [SfV 

mf l) =m- k+ K^(y k -n k ) 

p(' +1) =p--4 +1) 5f 1) [4' +1) f 

v£' +1) = V k ~ + J (y t - h(x k )) (y k - h{x k )f 

xN(jc t |m<°, Pf)Ax k . 

and set V k = Vf \ « 4 = mf \ P k = Pf . 

Algorithm 2: The Variational Bayesian Adaptive Gaussian Fil- 
ter (VB-AGF) algorithm 
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